Water plays a very important role in daily life and many physical, chemical, and biological processes.^[1,2] Although being one of the most common substances on earth, liquid water has many unusual thermodynamic properties. For example, at low temperatures, its isothermal compressibility and isobaric heat capacity increase (rather than decrease) upon cooling.^[1] Clarification of the local structure of liquid water is the key of understanding the origin of these anomalies.^[3,4] Furthermore, water molecules interact with each other and form a hydrogen-bond (H-bond) network. Any changes of the position or orientation of individual water molecules influence their neighbors strongly, yielding a collective motion in the H-bond network. The dynamic inhomogeneity of this network is also directly related to the local structures of liquid water.^[5,6]

It has been proposed that liquid water contain two types of local structures, and accordingly molecules in liquid water can be divided into two types: the high-density-liquid (HDL) and low-density-liquid (LDL), as revealed from the equilibrium distribution of some specific order parameters.^[7–9] The HDL’s local structure is disorder. There are four neighbors in its first shell and its fifth neighbor inserts into the gap between the first shell and the second shell. However, the local structure of LDL is ice-like with the fifth neighbor in the second shell. The presence of such two types can explain well the aforementioned unusual thermodynamic properties of water, and is thus widely accepted.^[1,3] Here we perform all-atom simulations and identify a third type of water molecules, which is characterized by a new type of (real and/or inherent) local structure and, in particular, by an over-coordinated H-bond structure.

The three water models we used are TIP4P/2005,^[10] SPC/E,^[11] and TIP4P-Ew.^[12] All the simulations were performed in isothermal-isobaric (NPT) ensembles, with those of TIP4P/2005 and SPC/E models performed by using the molecular dynamics package LAMMPS^[13] and those of TIP4P-Ew by GROMACS.^[14] Long-range solvers were used to compute the long-range Coulombic interactions.^[15,16] The simulation conformation is called “real structure”. To reduce the influence of thermal fluctuations on the local structures in liquid water,^[17] we used the conjugate gradient algorithm to minimize the potential energy of the whole system. The molecules conformation obtained in such a way is called “inherent structure”.^[18] In recent studies, such minimization method has been successfully employed to identify two different local structures of HDL and LDL.^[7–9]

The main analyses presented in this article are based on the TIP4P/2005 model. The simulation of this model includes two cases, containing 216 and 4000 molecules, respectively. The temperature of the 216-molecule simulations ranges from 200 K to 300 K. The simulation time is 200 ns at low temperatures. We repeated the simulation seven times at 200 K and 210 K. The shortest run is 100 ns at 300 K. The 4000-molecule simulations were performed at 300 K for 10 ns. The simulations for the TIP4P-Ew and SPC/E models contain 1000 molecules and were performed at 300 K for 20 ns.

In the simulation, the hydrogen bond is defined based on the distance R_OO between the two oxygens in different molecules (R_OO < 0.35 nm) and the angle θ_HOO between the covalent OH and OO bonds (θ_HOO < 30°).^[19]

The local structure index (LSI) of molecule i, I(i), is defined as follows.^[7] Assume the distances (r_j’s) between molecule i and its neighbours j’s can be sorted as r₁ < r₂ < r₃ ⋯ < r_n(i) < 0.37 nm < r_{n(i) + 1}, where 0.37 nm sets the starting distance of the second shell; we then have

Let P(t) denote the vector of probabilities of all the states at time t. According to the Chapman–Kolmogorov equation, we can obtain equation P(t) = T(t)P(0), where the matrix T(t) is the transition probability matrix. The transition probability matrix element T_ij(t) represents the probability of molecules jumping from state i to state j in time duration t. Assume there are N₁ molecules on state i at time t₀, and N₂ molecules among them jumped to state j at time t₀+t. The ratio between N₂ and N₁ then gives the value of T_ij(t). Also, the time evolution of P(t) is governed by equation Ṗ(t) = KP(t), where K is the transition rate matrix. The solution of this equation is P(t) = e^KtP(0). Comparing this solution and the Chapman–Kolmogorov equation, we find that the transition probability matrix T and transition rate matrix K at time t are related by T = exp(tK).^[20] The diagonalizable T can be factorized as T = QΛQ⁻¹. Each column of Q is the eigenvector of T and Q⁻¹ is the inverse of Q. Λ is a diagonal matrix containing all the eigenvalues of T. The logarithm of T, log(T), is equal to Q log(Λ)Q⁻¹. Then we can obtain the transition rate matrix .

We perform all-atom molecular dynamics simulations for the three microscopic water models at temperatures ranging from supercooled to ambient. By employing multiple order parameters to characterize the (real and/or inherent) local structure of water, we identify the presence of three types of molecules in liquid water: in addition to the well-known HDL and LDL,^[21–41] there exists a third type, characterized by an ultra-high density and over-coordinated H-bonds. When the temperature rises, the fraction of this third type increases, from a very small value (< 1%) in the deep supercooled region (200 K) to about 7% at the room temperature (300 K). The ultra-high-density-liquid (UHDL) molecules are surrounded by HDL molecules and distribute discretely in space. They decrease the probability of transition of the HDL molecules to LDL molecules and increase the probability of the reverse one, thus leading to the enhancement of the stability of the HDL molecules.

3.1. A new local structure

We start with observing the distributions of r₅ in the inherent structure for different temperatures (Fig. 1(a)), different pressures (Fig. 1(b)), and different water models (Fig. 1(c)). All the distributions have two peaks. The position of the left peak is smaller than r₅ = 0.3 nm, and that of the right peak varies slowly with the temperature or pressure. It is worth noting that the distributions obtained from different water models are almost the same, which proves that our results are independent of models. As shown in Appendix A, the bimodal characteristic of the distribution of r₅ can also be observed in the real structure.

Fig. 1.

Figure Option ViewDownloadNew Window

Fig. 1. Distributions of r5 in the inherent structure for different temperatures, different pressures, and different water models: (a) distributions of r5 in the inherent structure at different temperatures, calculated from the 216-molecule simulations at 1 atm; (b) distributions of r5 at different pressures, obtained by the 4000-molecule simulations at 300 K; (c) distributions of r5 for different models, given by the 1000-molecule simulations at 300 K and 1 atm.

To characterize the local structures in liquid water clearly, we calculate the joint distribution of r₅ and LSI in the inherent structure. As shown in Fig. 2(a), for water at 1 atm and 300 K, there exist three peaks, indicating the presence of three types of local structures of liquid water. The upper right peak, with large LSI and large r₅, corresponds to the LDL structure.^[7] Among the two peaks with small LSI, the upper left one, with larger r₅ and much larger height, consists the main part of the low-LSI component,^[7] and should thus correspond to the HDL structure; the lower left peak, with the smallest height, has the smallest r₅, representing a local structure with density even higher than that of HDL, i.e., the ultra-high-density liquid (UHDL). Note that the bimodal characteristic of the UHDL peak is caused by the fluctuation of the sixth neighbor around the cutoff at 0.37 nm. The projection of the joint distribution on the vertical axis goes back to the result in Fig. 1.

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. Identification and characteristics of the third type water molecules. (a) Joint I–r5 distribution at 1 atm and 300 K in the inherent structure, where the values of the probability density are given by the color bar on the right. The three peaks corresponding to LDL, HDL, and UHDL have been labeled with L, H, and U, respectively. The boundary between the UHDL and HDL is r5 = 0.305 nm and the boundary between the HDL and LDL is I = 1.35 × 10−3 nm2. (b) Average number of acceptor bonds of UHDL at various temperatures. In the inherent structure, the acceptor bond number is always close to three at all temperatures; in the real structure, due to thermal fluctuation, the acceptor bond number decreases when the temperature increases, but still remains bigger than two even at 300 K. (c) Schematic illustration of three species in the inherent structure, where the fifth neighbor of the LDL molecule locates beyond the 0.37 nm circle and that of the HDL molecule is in the gap between the first and second shells, while the UHDL molecule has five neighbors in the first shell.

We proceed to investigate the H-bond structure of UHDL molecules. Figure 2(b) shows the average number of the acceptor bonds of UHDL at different temperatures: in the inherent structure, the average acceptor bond number is always 3 at all temperatures; while in the real structure, due to the thermal fluctuation, the average acceptor bond number of the corresponding UHDL molecules decreases gradually, but remains much bigger than two even at 300 K. {Normally, water molecules have only two acceptors, but UHDL molecules have more than two. That is why we say the UHDL molecules are over-coordinated.

A schematic illustration showing the differences between the LDL, HDL, and UHDL is given in Fig. 2(c) (i) LDL is ice-like, with the fifth neighbor locating beyond the 0.37 nm cutoff circle; (ii) HDL has four neighbors in the first shell, and the fifth neighbor is in the gap between the first and second shells; (iii) UHDL, however, has five neighbors in the first shell, with no molecule in the gap. The LDL and HDL structures are consistent with the descriptions in Ref. [1]. Note that UHDL is different from the very-high-density amorphous ice (VHDA), which has four neighbors in the first shell and two interstitial molecules.^[42]

3.2. Roles of UHDL molecules

We label each molecule in the simulations according to its inherent structure. We next discuss the roles of the UHDL molecules. To this end, we first analyze the spatial correlations between the LDL, HDL, and UHDL molecules. Figure 3(a) gives a snapshot of liquid water at 300 K and 1 atm, where the green, white, and red spheres represent the oxygen atoms of LDL, HDL, and UHDL molecules, respectively. As shown in Fig. 3(a), the LDL and HDL molecules prefer to cluster with the same species, but the UHDL molecules are dispersed among the HDL. In order to confirm this observation, we calculate the percentage of species β in the first-shell neighbors of specie-α molecules, C_αβ, where α and β represent LDL, HDL, or UHDL (abbreviated as L, H, and U, respectively), and the concentration of species α in the whole system, C_α. By definition, we have ∑_α C_α = 1 and ∑_β C_αβ = 1. If C_αβ is larger than C_β, species α and β attract each other; otherwise,α and β repel each other. Figures 3(b) and 3(c) give, respectively, the values of the three C_α and the three C_αU at 1 atm and various temperatures between 200 K and 300 K, from which we can clearly see that HDL molecules attract UHDL molecules while the LDL molecules repel UHDL (because C_HU > C_U > C_LU). The C_UU is smaller than C_U, indicating that the UHDL molecules do not tend to form clusters. These two results together suggest that UHDL molecules probably serve as the clustering nuclei of HDL molecules.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. Concentrations (Cα) of and correlations (Cαβ) between different molecule types. (a) Simulation snapshot of 216 liquid-water molecules at 1 atm and 300 K, where the green, white, and red spheres represent the LDL, HDL, and UHDL molecules, respectively. These molecules are labeled according to their inherent local structures. Panels (b) and (c) show, respectively, the temperature dependence of Cα and at 1 atm, as given by the 216-molecule simulations.

The existence of UHDL also influences the transitions between the HDL and LDL molecules. To quantify such influences, we compare the transition matrix elements, T_HL (representing the probability of transition from the HDL to LDL) and T_LH (from LDL to HDL), of the molecules that are in the first shell of UHDL with those of the molecules that are not. As shown in Figs. 4(a) and 4(b), the HDL molecules in the first shell of the UHDL have a lower probability of jumping to LDL, and the first-shell LDL molecules have a higher probability of jumping to HDL, than their non-first-shell counterparts. In another word, the existence of UHDL makes its surrounding HDL more stable and LDL less stable.

Fig. 4.

Figure Option ViewDownloadNew Window

Fig. 4. Transition matrix elements. Panels (a) and (b) show, respectively, the TLH (the transition probability from L to H) and THL (the transition probability from H to L) at 1 atm and 300 K, given by the 4000-molecule simulations, where the open squares (circles) represent the HDL (LDL) molecules in the first shell of UHDL (HDL1 and LDL1) and the solid squares (circles) represent the non-first-shell HDL (LDL) molecules (HDL2 and LDL2).

All the transition probabilities and the corresponding transition rates between the different local structures are given in Fig. 5. The transition rate represents how fast one water molecule jumps from one state to another state. If the transition rate from a state to another is much faster than that between other states (orders of magnitude), then the stability of this metastable state is worse than that of the other states, so this metastable state is generally not treated the same as other states. As shown in Fig. 5(b), the transition rate from UHDL to HDL is of the same order as that from LDL to HDL, indicating that UHDL is a stable local structure.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (a) Transition probabilities (obtained by the 4000-molecule simulations at 300 K and 1 atm) and (b) transition rates between the different local structures, where the numbers in the circles in (b) represent the concentrations of the corresponding species. The convergent transition rates are obtained from the transition probability matrix at t = 0.8 ps.

We predicted the existence of a third type of molecules (UHDL) in liquid water by investigating its inherent structures, and its main characteristic remains in the real structure. One distinctive feature of such molecules is their over-coordinated H-bond structure, i.e., the existence of three acceptor bonds. The UHDL molecules disperse in space and are usually surrounded by HDL molecules. The existence of UHDL makes the HDL in contact more stable and LDL in contact less stable. The heterogeneity of liquid water is therefore greatly enhanced. Moreover, although the fraction of the UHDL molecules itself may be small, the fraction of their HDL/LDL neighbors can reach ∼ 35% (the proportion of UHDL molecules is ∼ 7% and each UHDL molecule has five neighbors, so the proportion of the neighbors of UHDL is ∼ 7% × 5) under ambient conditions. In this sense, the presence of UHDL has bigger influence than what its small fraction seemingly means. Its potential influence on the macroscopic properties of water will be explored in near future.